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THREE-DIMENSIONAL  ELECTROMAGNETIC  SCATTERING 
FROM  ARBITRARY  INHOMOGENEOUS  OBJECTS 


Introduction 

The  purpose  of  this  paper  is  to  describe  the  work  we  have  done  towards  de¬ 
veloping  a  versatile  and  efficient  method  for  calculating  electromagnetic  energy 
deposition  in  arbitrary  three-dimensional  inhomogeneous  objects. 

Given  an  incident  field  Emp  generated  from,  for  example,  an  antenna,  the 
problem  is  to  estimate  the  amount  of  electromagnetic  energy  deposition  in  a 
nearby  human  object.  While  this  is  a  classical  problem,  a  versatile  and  efficient 
method  for  solving  it  in  realistic  settings  is  not  generally  available  [10].  Some 
recent  research  in  this  area  can  be  found  in  [24,  25,  26]  and  their  references. 
Understanding  the  details  of  electromagnetic  deposition  in  humans  is  essential 
for  health  and  safety  exposure  considerations.  Knowledge  of  deposition  in 
humans  and  non-human  animals  is  very  important  to  medical  research  on  the 
bioeffects  of  radiation  exposure. 


Derivation  of  Model  Equations 

The  stated  problem  involves  solving  the  symmetric  Maxwell’s  equations: 

VxE  =  -  VJf 

at 

„  -  di)  -true 

V  x  H  =  — —  +  J* 
at 

V  D  =  ptrue 
V-B  =  ftf* 

for  the  electric  and  magnetic  field  intensities  E  and  H,  and  the  electric  and 
magnetic  flux  densities  D  and  B  everywhere  in  R3  including  the  object.  We 
assume  the  body  (scatterer)  occupies  a  bounded  region  V  in  free  space.  Here 
the  electric  sources  j?rue  and  Jtrue  are  respectively  the  charge  and  current 
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densities.  For  the  sake  of  symmetry,  the  fictitious  and  identically  vanishing 
magnetic  sources  p^ct  and  JAct  are  introduced.  All  dependent  (“tilded”)  vari¬ 
ables  are  assumed  to  be  functions  of  space  and  time. 

We  assume  the  object  is  linear,  isotropic  and  nondispersive.  Then 

D  =  eE 

B  =  /xH 

jtrue  =  aE 

where  e,  p,  and  a,  all  possibly  spatially-dependent,  are  respectively  the  permit¬ 
tivity,  the  permeability,  and  the  conductivity  of  the  body.  Outside  the  body, 
the  first  two  of  these  equations  hold  with  the  constant  free-space  parameters 
£o  and  fiQ  respectively. 

If  we  let  n( r)  =  p0  +  and  e(r)  =  e0  +  6e( r),  where  r  =  ( x ,  y,  z ),  then 
the  Maxwell’s  equations  can  be  written  as 

V  x  E 

V  x  H 

V-E 
V  H 

where 


J  = 

P  = 

Pm 

[f  we  further  assume  the  sources  vary  sinusoidally  with  time,  then  the 
problem  reduces  to  solving  the  time-harmonic  equations: 


V  x  E  = 

-jpup0H  -Jm 

(1) 

V  x  H  = 

jf3uj£0E  +  J 

(2) 

V-E  = 

p/£o 

(3) 

V  H  = 

Pm/ Po 

(4) 

P°  m  Jm 

dE  - 

e”ar  +  J 

P/So 
Pm /  Po 


c  .  t  true 

ptrue  -  V  •  8eE 

Pi T  -  V  •  6 pH 
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subject  to  appropriate  boundary  conditions  on  the  E  and  H  fields.  Here  the 
phasor  E  is  related  to  E  by  the  equation 

E(x,  y,  z ,  t)  =  5J[E(z,  y,  z)ejpojt} 

where  ft  =  ±1.  Similarly,  the  same  can  be  said  for  the  other  dependent 
variables.  Clearly,  we  must  have 


m  =  jpUlSpK  + 

(5) 

J  =  jpuSeE  +  Jtrue 

(6) 

p  =  Ptrue-V-Se  E 

(7) 

'm  =  PfJf  -  V  •  6 pH 

(8) 

If  the  electric  and  magnetic  intensities  E  and  H  satisfy  the  time  harmonic 
equations,  Equations  (l)-(4),  then 


V  J  =  -j/3ujp  (9) 

V  •  Jm  =  -jpUpm  (10) 

Or,  equivalently,  in  terms  of  the  original  variables, 

V  •  Jtrue  =  -j/3uptrue  (11) 

V-J£“  =  -jPvpg*  (12) 


These  are  the  equations  of  continuity.  Conversely,  if  E  and  H  satisfy 
Equations  (1)  and  (2),  then  Equations  (3)  and  (4)  are  automatically  satisfied, 
provided  the  equation  of  continuity  for  the  electric  sources  3true  and  ptrue  holds. 
Since  we  will  always  assume  this  to  be  true,  it  suffices  to  solve  only  Equations 
(1)  and  (2). 

Taking  the  “curl”  of  Equation  (1)  and  then  substituting  Equation  (2)  in 
the  result,  we  obtain  the  equation 

V  x  V  x  E  -  uj2£op0E  =  -j/3up0J  -  V  x  JTO  (13) 

The  formulation  up  to  this  point  has  been  quite  general.  However,  we 
will  now  make  the  simplifying  assumptions  that  p  =  pa  everywhere  inside  the 
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body.  It  follows  that  Jm  =  0  and  pm  =  0  everywhere,  as  J^cf  and  p^f1  vanish 
identically  in  reality.  Thus,  Equation  (13)  simplifies  to 

V  x  V  x  E  -  k20E  =  -jPup0J  (14) 


where 

K  =  W  £oHo 

If  we  let  G(r,  r')  denote  the  free-space  Dyadic  Green’s  Function,  so  that 

G(r,r')  =  (I  +  pVV)S(r,r') 


where 


r')  = 


g-jjSMr— r'| 

47t  |  r  —  r'  | 


is  the  Green’s  function  for  the  three-dimensional  scalar  Helmholtz  equation, 
then,  treating  the  right-hand-side  of  Equation  (14)  as  a  source  (even  though  it 
contains  the  unknown  field  E  through  the  definition  of  J),  it  can  be  shown  that 
the  electric  field  E  formally  satisfies  the  following  Fredholm  integral  equation 
(of  the  second  kind): 


E(r)  =  -jpujfio  [  G(r,r')  •  J(r')  dr' 

JVUVf 

=  -j(3 up0  [  G(r,  r')  •  [j(3u8e{ V)  E(r')  +  Jtrue(r')  ]  dr' 

Jvuvf 

=  -jfjun,  /  G(r,  r')  ■  \j0u6eV)  E(r')  +  J,r“(r')  ]  *' 

J  v 

-jPupof  G(r,  r')  •  3tTue(r')  dr' 

Jv< 

where  Vf  denotes  that  region  in  space  outside  the  body  V  where  the  “source” 
J  does  not  vanish.  Here,  we  have  used  the  fact  that  8e  =  0  in  Vc,  the  region 
external  to  V.  Since  we  are  assuming  the  body  medium  is  linear,  so  that 
Jtrue(r)  =  cr(r)  E(r),  it  follows  that 

E(r)  =  G(t,r>)-  [r(r')  E(r')  ]  *' 

J  V 

-jfaHof  G(r,r')-J‘r”(r')*'  (15) 

JVS 
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where 


/  x  c  f  v  j/Mr) 

r(r)  —  <?e(r) - — 

UJ 

Since  the  electric  field  E  is  identical  to  the  incident  field  Emc  when  there  is  no 
scatterer  (Se  =  0  and  a  =  0),  the  last  integral  in  Equation  (15)  must  be  the 
incident  field: 

Einc(r)  =  -jfano  [  G{ r,  r')  •  3true(r')  dr ' 

Jv§ 

Thus,  the  model  (integral)  equation  we  need  to  solve  is: 

E(r)  +  0O  f  G( r,  r')  •  [r(r')  E(r')  ]  dr1  =  Einc(r)  (16) 

J  v 

where 

Po  =  -u2Ho 

Singularity  of  the  Dyadic  Green’s  Function 

At  each  field  point  r  inside  the  region  V,  the  model  equation,  Equation 
(16),  becomes  singular,  since  G(r,r')  is  undefined  when  the  source  point  r1 
coincides  with  the  field  point  r. 

The  singularity  of  the  Dyadic  Green’s  Function  has  been  well-studied  ([17] ,  [18] , 
[15],  [20]).  The  singular  integral  in  Equation  (16)  can  be  handled  “as  usual” 
as  follows. 

(  Cfr.r'j  ■  [r(r')  E(r') ;  dr'  =  , J.im„  ^  G(r,r')  ■  (r(r')E(r')]  dr' + 

/  G(r,r')-[r(r')E(r')]rfr')  (17) 

where  Va  are  ever-shrinking  volumes  surrounding  the  field  point  r.  However, 
what  is  unusual  about  this  approach  is  the  fact  that  the  limiting  process  is  not 
uniform.  In  particular,  the  limiting  values  of  the  two  integrals  on  the  right 
will  generally  be  different,  depending  on  the  shapes  of  Vs.  Nevertheless,  the 
sum  of  the  two  will  always  be  the  same. 

If  the  Vs  are  chosen  to  be  spheres,  then  it  can  be  shown  ([17],  [15], [20],  [7]) 
that 

i*'  =  -hhfh  (is. 
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Thus,  Equation  (16)  becomes 

<*o(r)E(r)  +  P0  Upv(r)  =  E™(r)  (19) 

where 

«o(r)  =  1  +  ^  (20) 

and 

Upv(r)  =  PV  f  G(r,  r')  •  [r(r')  E(r')  ]  dr'  (21) 

J  V" 

Here  PV  stands  for  principal  value  integration  (using  shrinking  spheres).  This 
formulation  has  been  employed  by  Livesay  and  Chen  [12]  and  many  other 
researchers  ([5],  [18],  [4]). 

The  principal  value  integration  in  Equation  (19)  generally  cannot  be  eval¬ 
uated,  analytically  or  numerically.  An  equivalent  but  more  practical  approach 
is  to  re-formulate  the  singular  integral  in  Equation  (17).  Defining 

F(r)  =  r(r)  E(r), 


and  using  a  representation  of  E(r)  that  involves  vector  and  scalar  potentials, 
one  can  replace  this  singular  integral  by 

Ufn(r)  =  (I  +  iw)  •  Jv  g(r,  r')  F(r')  dr' 

Moreover,  Ufn(r)  can  be  expressed,  without  principal  value  integration,  as 
([8,  11,  19,  20]) 

Ufa(r)  =  I1(r)  +  I2(r)  +  I3(r)  (22) 


where 


I,(r)  =  f  G(r,r')-F(r')*' 

J  V  —  Vrp 

I2(r)  =  /  G(r,r').F(r')-G„(r,r')-F(r)*' 

JVt 


J_  f  n(rQ  (r  -  rQ 

kl  J  St  47t  |  r  -  r'  |3 


F(r)  dS' 
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Here,  the  region  Vt  is  any  finite  (not  necessarily  infinitesimal)  region  enclosing 
r  and  is  contained  in  V.  St  is  the  surface  area  of  Vt-  The  outward  normal  of 
St  at  a  point  r'  on  St  is  denoted  by  n(r').  Finally,  G0  is  the  so-called  static 
Dyadic  Green’s  Function  defined  by 


G„(r,r')  =  iw9o(r,r') 


where 


9o{  r,  r')  = 


47 r 


1 


r  —  r' 


In  their  derivation  of  the  representation  above,  Fikioris  [8],  Lee  [11]  ,  and  Wang 
([19],  [20])  require  F  to  satisfy  a  Holder  condition.  F  satisfies  a  Holder  condi¬ 
tion  at  a  point  r  if  there  are  three  positive  constants  A,  B,  and  C  such  that 


|  F(r)  -  F(r')  |<  A  |  r  -  r'  \B 


for  all  r'  satisfying  |  r  —  r'  |  <  C 
It  is  interesting  to  note  that 


PV  J  G0( r,  r')  •  F(r)  dr'  ±  ±  j  n(r')V'5o(r,  r')  dS'  •  F(r) 

_  J_  r  n(r')  (r  -  r') 
k%  J sT  -47T  |  r  —  r'  I3 


uO  JST 

=  W 


This  can  be  verified,  for  example,  by  assuming  Vt  is  a  sphere  centered  at  r. 
For  in  this  case, 

PV  /  G0(r,r')-F(r)  dr'  =  0 

J  Yx 

but 


hi  JsT  An  |  r  —  r' 


Consequently, 

Ufn(r)  ±  f  G(r, r')  •  F(r')  dr'  +  PV  /  G(r, r')  •  F(r')dr' 

J  V — Vj7  J  Yj1 

=  PV  J^G(r,r')  ■  F(r')  dr' 
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This  implies  that  the  relationship 


(I  +  ^VV)  ■  jv  g(r,  r')  F(r')  dr'  =  PV  J^l  +  ± VV)  •  g( r,  r')  F(r')  dr '  (23) 

is  at  best  formal.  The  left  side  of  Equation  (23)  is  the  correct  term  to  use, 
though  it  is  numerically  undesirable.  Its  equivalent  three-term  form  in  Equa¬ 
tion  (22),  while  more  complicated,  is  numerically  better  behaved.  Finally,  the 
right  side  of  Equation  (23)  is  incompletely  defined,  requiring  additional  qual¬ 
ifications  (principal  volume  being  used  and  its  associated  correction  term)  for 
its  correctness. 

One  can  also  show  in  the  special  case  in  which  VT  is  a  sphere  that  I2  will 
approach  zero  as  the  radius  of  the  sphere  approaches  zero,  so  that  a  result 
implicit  in  Equation  (19)  is  recovered: 

(I  +  pVV)  •  Jv  g( r,  r')  F(r')  dr'  =  PV  ^  G(r,  r')  •  F(r')  dr'  -  j^F(r)  (24) 
Or,  equivalently, 

Ufn(r)  =  Upv(r)  -  gijF(r) 

In  short,  using  the  representation  of  the  singular  integral  given  in  Equation 
(22),  we  obtain  an  alternate  representation  of  the  model  equation  given  in  (19) 
but  which  does  not  require  making  a  principal  volume  integration: 

E(r)  +  /?0Ufn(r)  =  Einc(r)  (25) 

This  has  exactly  the  same  form  as  Equation  (16)  except  the  integral  is  now 
unambiguously  defined. 


Solution  Method 


To  solve  the  integral  equation  (19)  or  its  alternate  form  in  equation  (25), 
we  use  the  classical  Moment  Method  (MM)  [9].  In  this  report,  we  will  concen¬ 
trate  on  solving  equation  (19).  Because  we  do  not  assume  the  scatterer  to  be 
homogeneous,  E(r)  is  generally  not  divergence-free.  (See  Equations  (3)  and 
(7).  )  As  we  will  see  in  the  next  section,  it  is  consistent  with  our  approach 
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if  the  field  being  sought  is  divergence-free.  Thus,  instead  of  applying  MM  to 
equation  (19),  we  will  apply  it  to  the  following  equivalent  integral  equation: 

«i(r)E(r)  +  Po  PV  f  G(r,  r')  •  [n(r')  E(r')]  dr’  =  Einc(r)  (26) 

J  v 

where 


ax(r) 

di(r) 

n(r) 

E(r) 


=  a0(r)di(r) 

£o 

e0  +  r(r) 

=  di(r)r(r) 
=  E(r)/«i(r) 


Clearly,  di(r)E(r)  =  E(r)  and  rx(r)E(r)  =  r(r)E(r)  =  F(r).  More  impor¬ 
tantly,  it  can  readily  be  verified  that  E  is  now  divergence- free  (V  •  E  =  0.) 

In  MM,  the  unknown  field  E  is  assumed  to  be  expandable  in  a  family  of 
basis  functions  {fn}: 


E(r)  =  £>nfn(r)  (27) 

Tl~  1 

where  {an}%=1  are  complex  constants  to  be  determined.  When  Equation  (27) 
is  substituted  into  Equation  (26),  we  readily  obtain 

•£on{tti(r)iVi(r)  +  /?0Un(r)}  =  Einc(r)  (28) 


U„(r)  =  PV/  G(r,r,)-[r1(r')^(r')]dr'  (29) 

Starting  with  Equation  (28),  there  are  two  common  methods  to  determine 
the  N  complex  constants  {an}^=1.  One  can  evaluate  Equation  (28)  at  N 
distinct  points,  r  =  r i,i  =  1, ...,  N,  to  obtain  a  system  of  N  linear  equations 
in  N  unknowns.  This  is  the  so-called  point-matching  method.  Alternately, 
we  can  integrate  Equation  (28)  with  each  of  the  N  basis  functions  fn  in  turn 
to  get  again  a  system  of  N  linear  equations  in  N  unknowns.  This  is  the  so- 
called  Galerkin’s  Method.  This  is  the  method  we  have  used  in  this  study.  The 
optimality  of  the  Galerkin’s  Method  has  recently  been  discussed  in  [21]. 
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Two  issues  remains  to  be  addressed  to  complete  the  description  of  the 
solution  method.  One  is  the  choice  of  basis  functions  {fn}  and  the  other  is  a 
method  to  evaluate  the  integral  defining  Un(r)  in  Equation  (29). 


Edge-based  vector  basis  functions 

All  numerical  solutions  are  approximations  of  the  exact  solutions.  Thus, 
even  though  exact  solutions  to  Equations  (1,  2)  will  automatically  satisfy 
Equations  (3,  4)  in  light  of  the  Equations  of  Continuity,  Equations  (11,12), 
numerical  solutions  to  Equations  (1,  2)  may  be  far  from  satisfying  Equation 
(ii). 

The  requirement  that  Equation  (3)  be  satisfied  can  readily  be  shown  to  be 
equivalent  to  the  condition  V  •  E  =  0.  Since 

N 

E(r)  =  Y,  anfn(r), 

n— 1 

this  condition  is  automatically  satisfied  if  V  •  fn  =  0  for  each  f„.  In  finite 
element  methods,  a  well-known  class  of  functions  with  this  property  is  the 
class  of  Whitney  functions  of  degree  1  ([23],  [3]).  We  adopt  this  family  for 
our  use  here.  For  sake  of  completeness,  we  will  define  this  family  and  mention 
some  of  its  properties  relevant  to  our  numerical  procedure. 

In  the  numerical  solution  of  Equation  (19),  we  approximate  the  region  V 
occupied  by  the  body  by  a  family  of  Nt  disjoint  tetrahedra  {T*}^.  Since  there 
are  six  edges  to  a  tetrahedron,  there  will  be  a  large  number,  say  Ne,  of  edges 
in  the  entire  system.  (There  is  generally  no  formula  that  relates  Nt  to  Ne. 
Two  tetrahedra  can  have  12,  11,  or  9  edges  depending  on  whether  they  have 
no,  1  or  three  edges  in  common.  In  our  discretization,  we  do  not  allow  two 
tetrahedra  to  have  touching  faces  without  them  being  the  same.)  Similarly, 
since  there  are  four  vertices  (nodes)  to  a  tetrahedron,  there  will  also  be  a  large 
number,  say  Nn,  of  nodes  in  the  system. 

To  each  edge  in  the  system  is  associated  a  unique  Whitney  function.  If 
e  =  (vj,  V2)  denotes  an  edge  in  the  system  joining  the  node  Vj  to  the  node 
V2,  then  the  vector- valued  Whitney  function  We  associated  to  the  edge  e  = 
(vi,  V2)  is  defined  by 

W,(r)  =  AVl(r)VAV2(r)  -  AV2(r)VAVl(r) 
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where  Av,  (r)  is  the  barycentric  function  associated  to  the  node  v,,  (i  =  1,2), 
and  is  the  “simplest”  piecewise  linear  function  such  that  AVi(vj)  =  b{ ,  j  = 
1, . . . ,  Nn.  In  particular,  if  Vj,  (i  —  1, . . . ,  4)  represent  the  four  vertices  of  a 
tetrahedron  T,  then  AVi  has  the  following  simple  representation  in  T: 


Aj(r) 


—  G7; 


— 1 


d=  bj  ■  r  +  aj 


(30) 


where 


and 


Vi  v2  v3  v4 
1111 


<*  =  (4,44,4] 

Furthermore,  if  e  =  (vi,  v2)  in  an  edge  of  T,  a  representation  of  We  in  T  is  of 
the  form  (using  notations  in  Equation  (30)): 

Wer(r)  =  aj2  +  bf,2  x  r  (31) 


where 


and 


ai'2  =  afb2  -  o2  b( 


b^2  =  bf  x  b^ 

With  some  algebraic  manipulations,  one  can  also  show  that 

v3  x  v4 


al,2 
/  _ 


bl, 


det(X) 

V4  ~  V3 

det(X) 


where  det(X)  =  (v2  -  v4)  •  [(v3  -  Vi)  x  (v4  -  Vi)]  is  the  determinant  of  X. 
Hence,  the  Whitney  function  associated  with  the  edge  e  =  (vi,  v2)  has  a  rather 
simple  representation  in  the  tetrahedron  T  whose  vertices  are  vi,v2,v3,  and 
v4: 


wHr)  =  ^(V3xv4  +  (v4-v3)xr) 

(v3  -  r)  x  (v4  -  r) 
det(X) 


(32) 

(33) 
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It  should  be  noted  that  while  a^2,  b^2  and  therefore  Wj  depend  on  all  four 
vertices  Vi, . . . ,  v4  of  T,  they  are  invariant,  as  they  should,  under  the  ordering 
of  v3  and  v4. 

The  following  properties  of  We  are  especially  relevant  to  our  numerical 
method.  Firstly,  the  We  is  divergence-free,  i.e.,  V  •  Wj  =  0  for  any  r eT.  This 
follows  directly  from  Equation  (31). 

Secondly,  it  can  be  shown  algebraically  that  the  tangential  component  of 
Wj  on  each  of  the  two  faces  of  T  containing  the  edge  e  is  dependent  only 
on  the  vertices  making  up  that  face,  and,  on  the  other  two  faces,  the  tangent 
component  of  Wj  is  identically  zero.  Since  We  vanishes  identically  on  any 
tetrahedron  not  containing  the  edge  e,  the  tangential  component  of  We  is 
continuous  across  all  faces  of  the  tetrahedron. 


Evaluation  of  Un 

To  evaluate  the  integral  in  Equation  (29),  we  re-write  it  as  follows: 

U„(r)  =  PV^G(r,r')[r1(r')f„(r')]<ir' 

=  IJ(r)  +  IJ(r) 


where 

JV-VT(  r) 

I£(r)  =  PV/  G(r,  r')  •  7i(r,)fn(r/)  dr' 

Here  Vr(r)  is  the  unique  tetrahedron  containing  r  in  its  interior.  (In  the 
Galerkin’s  Method,  Un(r)  is  eventually  multiplied  (dot-producted)  with  the 
basis  functions  (Whitney’s  functions)  and  integrated  numerically  over  tetrahe- 
dra.  The  second-order  numerical  integration  method  employed  only  requires 
the  evaluation  of  the  dot-product  at  interior  points  of  tetrahedra.) 

The  integral  in  I2  is  generally  difficult  to  evaluate  over  tetrahedra.  For  this 
reason,  we  make  the  simplifying  assumption  that 

I£(r)  «  PV  f  G(r,r')-ri(r')fn(r')  dr' 

JB(r,req) 
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where  B( r,  req)  is  a  sphere  centered  at  r  with  radius  req  given  by 

=  (£|V4.<r)|)* 

With  this  radius,  the  volume  of  the  sphere  is  exactly  the  same  as  the  volume 
of  Vt(t),  denoted  here  by  ||Vr(r)||.  Thus 

Un(r)  «  I?(r)  +  PV  f  G( r,  r')  •  n(r')  fn(r')  dr'  (34) 

JB(r,req) 

In  many  problems,  we  can  reasonably  assume  the  region  V  is  piecewise 
inhomogeneous.  In  these  cases,  we  can  clearly  ensure  each  tetrahedron  is 
within  a  region  of  constant  T\.  Since  the  basis  functions  fn(r')  in  this  study 
are  of  the  form  an  +  bn  x  r',  the  second  integral  in  Equation  (34)  can  be 
evaluated  explicitly: 


where 


PV  [  G{ r,  r')  •  ri(r')  fn(r')  dr' 

JB(r,req) 

Ti(r)  PV  [  G( r,  r7)  •  (a„  +  bn  x  r')  dr' 

JB(r,req) 


2T(reg)ri(r) 
3  kl 


(an  +  bn  x  r) 


T  (r^)  =  (1  +  jPk0req)e~jl3koreq  —  1 
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Results 


To  test  our  approach,  we  calculated  the  E  field  in  a  homogeneous  sphere 
irradiated  by  a  plane  wave  and  compared  the  result  with  the  analytic  solution 
(Mie).  We  did  this  for  several  different  sets  of  parameter  values  and  discretiza¬ 
tions. 


Test  Group  I. 

Our  first  sets  of  parameters  are  those  reported  in  [15].  In  particular,  we 
studied  two  different  problems:  one  is  for  a  transparent  sphere  in  which  the 
conductivity  a  is  zero  and  the  relative  dielectric  constant  e(r)/e0  is  identically 
one,  and  the  other  is  for  a  translucent  sphere  in  which  a  =  0.015625  mhos 
per  meter  and  the  relative  dielectric  constant  e(r)/ec  is  identically  1.015625. 
In  both  cases  we  assume  the  incident  field  is  a  plane  wave  propagating  along 
the  z-axis  and  polarized  in  the  x-direction.  The  frequency  is  taken  to  be  1000 
megahertz  and  the  amplitude  is  1  volt  per  meter.  To  keep  the  problem  small, 
we  assume  the  sphere  has  a  5  cm.  radius. 

For  the  assumed  frequency,  the  free  space  wavelength  is  approximately 
0.3  (m).  If  we  use  the  commonly  quoted  rule  of  thumb  of  20  points  per 
wavelength,  the  mesh  size  is  approximately  0.015  (m)  in  each  of  the  three 
directions.  A  tetrahedron  with  three  sides  parallel  to  the  three  rectangular 
coordinate  axes  and  each  with  length  of  0.015  (m)  has  a  volume  of  5.625  x  10-7 
(m3),  requiring  approximately  930  tetrahedra  to  fill  out  the  given  sphere.  If  we 
use  10  points  per  wavelength  instead,  then  the  volume  of  a  typical  tetrahedron 
is  4.5  x  10-6  (m3),  requiring  approximately  116  tetrahedra  to  fill  out  the  sphere. 
To  keep  our  test  problem  small,  we  have  discretized  the  sphere  using  a  smaller 
number  of  tetrahedra,  which  in  essence  amounts  to  using  a  smaller  number  (10 
approximately)  of  points  per  wavelength.  We  have  designed  and  implemented 
two  different  discretization  methods  (D1  and  D2)  for  the  sphere.  Method  D1 
was  used  for  the  transparent  case  and  D2  was  used  for  the  translucent  case. 
The  major  characterisitics  of  the  grids  resulting  from  these  discretizations  are 
summarized  in  Table  1. 

The  discretizations  we  have  chosen  for  simplicity  are  relatively  coarse,  as 
one  can  see  from  the  following  volume  calculation.  The  volume  of  our  5  cm- 
radius  sphere  is  5.23  x  10-4  (m3),  but  the  total  volume  occupied  by  all  the 
tetrahedra  in  either  the  transparent  or  the  translucent  case  is  approximately 
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4.0  x  10-4  (m3),  which  is  only  75%  of  the  actual  volume.  For  comparison,  the 
fourth  column  (“Preferred”)  in  Table  1  shows  an  example  of  a  discretization 
using  method  D2  in  which  97.59%  of  the  volume  of  the  sphere  was  achieved 
by  all  the  tetrahedra. 

Figures  1-9  show  graphically  the  results  of  these  discretizations.  As  one 
may  see  from  the  plots  of  the  shrunken  tetrahedra,  the  main  difference  between 
discretizations  D1  and  D2  is  that  the  tetrahedra  in  D2,  especially  those  near 
the  boundary,  are  more  regular  (not  as  flat)  than  those  in  Dl.  Incidentally,  the 
volumes  of  the  tetrahedra  in  the  translucent  case  have  a  narrower  distribution 
than  that  in  the  transparent  case.  However,  this  is  not  due  to  the  difference 
in  the  two  methods,  but  rather  to  the  choice  of  parameters  used  to  generate 
these  tetrahedra. 

The  “Preferred”  case  (Figures  7-9)  shows  in  particular  that  the  discretiza¬ 
tion  method  D2  will,  in  the  limit,  reproduce  the  sphere. 

Figures  10-11  show  the  magnitude  of  the  x-component  of  the  total  field 
inside  the  sphere  along  the  z-axis.  They  agree  quite  well  with  the  theoretical 
results,  despite  the  somewhat  coarse  discretizations. 

Test  Group  II. 

To  test  the  convergence  of  our  method,  we  refined  the  discretization  for  the 
translucent  case  further.  Figure  12  shows  the  result  of  using  288  tetrahedra. 
This  result  was  unexpected,  as  it  is  clearly  not  as  good  as  that  of  using  120 
tetrahedra  (Figure  11).  We  also  tried  using  400  tetrahedra,  the  result  (not 
shown)  did  not  get  better.  We  suspect  this  non-convergence  may  have  to  do 
with  our  handling  of  the  singularity  of  the  dyadic  Green’s  function. 

We  also  tested  the  sensitivity  of  our  method  to  parameters  of  the  model. 
Figure  13  shows  the  result  of  using  a  relative  dielectric  constant  of  2.015625 
(one  unit  more  than  what  was  used  to  produce  Figure  11.)  While  this  result 
is  apparently  not  as  good  as  that  shown  in  Figure  11,  it  is  not  unreasonable, 
considering  the  gross  discretization  and  the  increased  discontinuity  of  dielectric 
constant  (between  free  space  and  the  scatterer.) 

Finally,  we  tested  one  case  in  which  the  frequency  was  increased  from  1 
GHz  to  10  GHz.  The  result  is  shown  in  Figure  14.  Here,  the  model  has  consis¬ 
tently  under-estimated  the  actual  value,  though  we  can  see  a  mild  convergence 
towards  the  real  solution  as  the  discretization  is  refined. 
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Discussion 


Our  preliminary  tests  clearly  indicate  that  we  should  next  address  the 
convergence  problem.  This  is  likely  to  involve  better  handling  of  the  singularity 
of  the  dyadic  Green’s  function,  as  this  is  the  principal  approximation  used  in 
this  method. 


TABLE  1.  RUN  PARAMETERS 

Parameters 

Transparent 

Translucent 

Preferred  Grid 

Discretization 

Method 

D1 

D2 

D2 

Sphere  Radius 
(m) 

0.05  (m) 

0.05  (m) 

0.05  (m) 

No.  Tetrahedra 

100 

120 

3648 

No.  Edges 

323 

400 

11,400 

No.  Nodes 

98 

103 

3,127 

Total  Volume 
(m3) 

4.00  x  10-4 

4^ 

o 

00 

X 

o 

1 

5.11  x  10"4 

Ave  Tetra.  Vol 
(m3) 

4.0  x  10“y 

3.4  x  10~e 

1.4  x  10~7 

Min.  Tetra  Vol 
(m3) 

1.5  x  hH 

1.9  x  10"e 

7.4  x  10-8 

Max.  Tetra  Vol 
(m3) 

9.0  x  10"® 

6.7  x  10~b 

4.0  x  IQ-7 

Conductivity, 
a  (mhos/m) 

0.0 

0.015625 

Rel.  Dielect, 

(  e/eQ  ) 

1.0 

1.015625 

Frequency 

(MHz) 

1,000 

1,000 

Incident  Field 

Ex  (volt/m) 
Ej,(volt/m) 
(volt/m) 

1.0 

0.0 

0.0 

1.0 

0.0 

0.0 
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Conclusion 


In  this  report  we  describe  the  progress  we  have  made  towards  developing 
a  3D  EM  interior  scattering  model  to  predict  energy  deposition  in  realistic 
biological  media.  The  volume  integral  approach  taken  is  a  natural  and  math¬ 
ematically  sound  method  to  solve  this  problem.  While  we  have  made  much 
progress,  there  are  obviously  several  areas  that  required  further  analysis.  These 
include 

1.  Refine  singular  integral  calculations. 

2.  Upgrade  our  current  method  or  develop  new  volume  integral  equation 
method  to  solve  our  problem. 

3.  Conduct  numerical  analysis  on  the  method. 

4.  Continue  validating  the  model  with  known  solutions  and  experimental 
data. 

Associated  with  these  analyses  is  the  important  question  of  computational 
efficiency.  Clearly,  it  is  paramount  that  our  model  be  accurate.  However,  for 
a  model  to  become  a  useful  tool,  one  must  be  able  to  apply  it  to  realistic 
situations.  In  our  case,  this  means  using  the  model  on  scatterers  of  reasonable 
size  (not  just  small  isolated  organs).  Therefore,  in  the  future,  we  need  to 
conduct  the  following  closely-related  mathematics/computer  sciences  research: 

1.  Develop  methods  to  handle  larger  systems. 

2.  Explore  efficient  ways  to  model  pulses. 

3.  Explore  the  possibility  of  using  parallel  computing  to  speed  up  the  cal¬ 
culations. 

Finally,  much  research  has  been  and  is  continued  to  be  done  to  solve  exte¬ 
rior  scattering  problems  due  to  its  obvious  military  significance.  The  interior 
scattering  problem  that  we  are  interested  in,  on  the  other  hand,  has  received 
relatively  little  attention.  Our  emphasis  here  is  to  develop  a  model  of  interior 
scattering  based  on  rigorous  mathematics  to  address  realistic  biological  prob¬ 
lems  (3-dimensional  inhomogeneous  scatterers) .  We  believe  we  are  heading  in 
the  right  direction. 
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Figure  1 .  Discretization  (D1)  of  a  Sphere 


Figure  3.  Tetrahedron  Volume  Distribution  (D1) 


Tetrahedron  Volume  1Q-6 
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Figure  5.  Discretization  (D2)  of  a  Sphere  (Shrunk) 


Figure  6.  Tetrahedron  Volume  Distribution  (D2) 
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Figure  7.  Preferred  Discretization  of  a  Sphere 
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Magnitude  of  x-component  of  E  field 


Figure  10.  Transparent  Case 
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Figure  12.  Translucent  Case 
Magnitude  of  Ex  along  z-axis 
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Figure  13.  Sensitivity  /  er 
Magnitude  of  Ex  along  z-axis 
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Figure  14.  Sensitivity  /  mhz 
Magnitude  of  Ex  along  z  —  axis 


-  =  true  solution 
*  =  8  tetra 
o  =  64  tetra 
x  =  120  tetra 


